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Abstract 

Within the coexistence region between liquid and vapor the equilibrium pressure of a simulated 
fluid exhibits characteristic jumps and plateaus when plotted as a function of density at constant 
temperature. These features exclusively pertain to a finite-size sample in a periodic box, as they 
are washed out in the bulk limit. Below the critical density, at each pressure jump the shape 
of the liquid drop undergoes a morphological transition, changing from spherical to cylindrical to 
slab-like as the density is increased. We formulate a simple theory of these shape transitions, which 
is adapted from a calculation originally developed by Binder and coworkers [J. Chem. Phys. 120 , 
5293 (2004)]. Our focus is on the pressure equation of state (rather than on the chemical potential, 
as in the original work) and includes an extension to elongated boxes. Predictions based on this 
theory well agree with extensive Monte Carlo data for the cut-and-shifted Lennard-Jones fluid. We 
further discuss on the thermodynamic stability of liquid drops with shapes other than the three 
mentioned above, like those found deep inside the liquid-vapor region in simulations starting from 
scratch. Our theory classifies these more elaborate shapes as metastable. 
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I. INTRODUCTION 


Statistical mechanics gives a general prescription for extracting equilibrium properties 
from a system Hamiltonian. To this aim, the system partition function should be determined 
in any of a number of equivalent descriptions (ensembles), a task requiring the evaluation 
of a complicate multidimensional integral. In practice, this can only be accomplished in 
a few (mainly trivial) cases, which thus obliges one to go over to partial calculations and 
approximate theories, often of limited scope. With the advent of computer simulation, the 
analytically intractable program of statistical mechanics could eventually be attacked and 
its solution hnally came within reach, at least for hnite, not too large systems. 

Occasionally, however, numerical simulation may lead to misconceptions. An example 
is the loop found in the pressure equation of state constructed by simulation within the 
coexistence region of liquid and vapor. Similar loops are found in the van der Waals theory 
of condensation, which uses the double-tangent construction to make the chemical potential 
everywhere concave as a function of pressure at constant temperature. As a result, the 
so-called metastable and unstable branches of the original equation of state are thrown out 
as unphysical. It is so customary to associate condensation with van der Waals theory, that 
it would tempting to interpret the loops canonical-ensemble simulations always produce in 
isotherms of intensive variables as van der Waals loops and, therefore, to read them as a 
sign of the entrance of vapor in a metastable regime. However, as originally remarked in 
Refs. [U |2], the non-concave regions observed in the equation of state of a hnite system 
are fully equilibrium features arising from the use of periodic boundary conditions in the 
simulation. In particular, the hrst inversion of concavity encountered corresponds to the 
hrst occurrence of liquid-vapor separation. We have recently showed that, despite their fake 
character, one can take advantage of these pressure loops to extract, by plain thermodynamic 
integration, the right liquid-vapor coexistence parameters Em¬ 
in a series of papers EEHZ], Binder and coworkers carried out extensive grand-canonical 
Monte Carlo (MC) simulations of the Lennard-Jones (LJ) huid in a periodic cubic box, 
biasing the sampling in such a way as to constrain the system to stay in the two-phase 
region. They found a whole sequence of so-called shape transitions between various “phases” 
differing in the shape of the liquid droplet coexisting with vapor. Specihcally, upon increasing 
the system density p the liquid drop changed from spherical (“sph”) to cylindrical (“cyl”) 
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to slab-like (“slab”). Upon increasing p further, the reversed sequence of transitions was 
observed, with interchanged roles between liquid and vapor. At each rearrangement of 
the liquid-vapor interface, the chemical potential undergoes a sharp drop, followed by a 
density interval where it stays nearly constant (this region will be referred to in the following 
as a “plateau”). MacDowell et al. have successfully analysed these shape transitions by 
means of a capillary-drop theory UM- In a further series of paper [HHIQ], Binder and his 
group exploited their extremely accurate simulation results to extract information about the 
interface free energy of curved liquid-vapor interfaces, and thus have access to the Tolman 
length nn, an all-useful parameter in nucleation theory [T2HI1] (the only caution here would 
be that in Refs. [HHIQ] the interfaces were taken at full equilibrium rather than under the 
metastable conditions typical of nucleation experiments). 


We hereby consider the phenomenon of condensation in canonical-ensemble simulations 
in the light of a yet different theory which, at variance of van der Waals theory, right from 
the outset takes into account the periodic repetition of the system in space. The present 
theory originates from the few-line calculation presented in Sect. II-A of Ref. [5], but now 
putting the emphasis on the pressure - rather than the chemical potential - equation of 
state, since it is the pressure that is directly accessible in a canonical-ensemble simulation. 
The theory is further refined by an extension to elongated-cubic simulation boxes. As a 
matter of principle, geometric shapes other than spherical, cylindrical, and tetragonal may 
also occur for the liquid drop in equilibrium with vapor, and this possibility will be looked 
closely within our theoretical framework. 


The remainder of the paper is structured as follows. In Sec. II, we expose the details 
of our theory. In Sec. Ill, we compare simulation results for the cut-and-shifted Lennard- 
Jones potential with theoretical predictions, discussing the relative stability of spherical, 
cylindrical, and slab-like drops as a function of the box aspect ratio. The stability of other 
interface shapes, like those observed in simulation runs performed independently near the 
cylinder-slab transition density |1|, is also studied. Finally, we give our conclusions in Sec. 
IV. 
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II. THEORY 


In order to determine the pressure behavior of a fluid near above the vapor coexistence 
density, the idea originally put forward in Ref. [S] is to compare the Helmholtz free energy F 
of the homogeneous vapor with that of various competing heterogeneous “phases”, differing 
in the shape of the liquid drop in thermal equilibrium with vapor. Since the thermodynam¬ 
ically stable phase for hxed temperature T, volume V, and particle number N is the one 
with minimum F, the drop adopts the shape which makes the total free energy as small as 
possible, consistent with the amount of liquid present at the given density p = N/V. 

If the focus is on the pressure equation of state, rather than on the chemical potential, 
it is convenient to express the free energy in terms of the specihc volume v = 1/p. This 
choice offers a number of practical advantages as it will be made clear below. Let us hrst 
estimate the free energy of the homogeneous vapor as a function of v at hxed T. Denoting 
/(T, v) = F/N the free energy per particle, a general relation valid up to second-order terms 
in the deviations AT = T — Tq and Av = v — Vq from a given state point Sq = (Tq, Vq) is: 

f = fo- ^AT - PAv - %AT^ - ^AvAT + -^Av'^ , (2.1) 

N 2T Kt 2vKt 

with /o = /(To, uo)- In fhe latter equation S is the entropy, cy is the const ant-volume specihc 

heat, ap is the isobaric expansion coefficient, and Kp is the isothermal compressibility, all 

computed at Sq. Choosing Sq to be the condensation point at temperature T, for AT = 0 

and Av = v — = 1/p^ being the vapor specihc volume at coexistence) it follows from 


Eq. (2.1) that 


AThom = F - F^ = Pv(nv - v)N + ^ {v - , (2.2) 

denoting Py and Ky the condensation pressure and isothermal compressibility of the bulk 
vapor at py (in the following, we use “v” and “ 1 ” subscripts to denote bulk properties of the 
coexisting vapor and liquid). 

For a vapor system in equilibrium with a spherical drop of liquid hosting Ni atoms, in 
the capillarity approximation the cost of droplet formation is 


ATsph = Tv+i -Fy = {N- N,)fy + Ai/i + {367ry/^'y{NrVif/^ - Nfy 

= iVi(/i - /v) + (367r)i/3^(iVini)2/3 , (2.3) 


where 7 is the surface tension at temperature T (he., the free energy of the planar liquid- 
vapor interface). In deriving Eq. (2.3), curvature corrections to surface tension [Hj as well as 
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anisotropy effects [I5l [16] have been neglected. Assnming for the vapor and liqnid fractions 
the same chemical potential and pressnre as in the bnlk limit, the difference /i — /v becomes 
— ni). Moreover Nv = NiVi + {N — N{)Vv, as no volnme is attached to the interface 
(this is to be contrasted with the lever-rnle estimate of Vi made in Ref. [5], which implicitly 
assnmed zero adsorption for the interface). In conclusion, we get: 


AFsph = Pv(^^v - v)N + (367r)^/'’7n( 


1 / 3 ^,^ 2/3 / 


Vi 


2/3 


/Y2/3 _ 


(2.4) 


From Eqs. (2.2) and (2.4|), we readily derive the pressures of the two “phases”: 
P-Pv 


.Phom Pv 


K.P 


and = (2.5) 

3 Pi - Pv V P - Pv / 


At variance with Ref. [S] Phom is a non-linear function of p, which makes it better suited to 
reproduce the true pressure behavior close to (see Fig. 1 below). By the way, the (p —Pv)^ 


term in Phom may not be the exact one since Eq. (2.2) is only a second-order truncated 
expansion. 

When the liquid drop is a cylinder extending along the shorter edge of a cuboidal box, 
hence of length = Ly = L^ja = {V/ayP (where a > 1 is the box aspect ratio), the 
free-energy excess over the vapor at p^ becomes: 


AF.,,1 = jV|(/i - /.) + 2ir7rc,|L., 


( 2 . 6 ) 


where the radius rcyi follows from 


Trr^yiLj, = N\V\ with W = 


Vv — V 


N. 


Uv - Vl 


Upon inserting Eq. (2.7) into (2.6), we obtain: 


AP„i = P.(n. - v)N + 27iP^'ya-P^p-P^ ( 

VPi -Pv/ 


thus yielding 


Pcyi = Pv — 


- Av / Vy-vA ^_i /3 

3(jl/6„5/6 _ y / 


(2.7) 


( 2 . 8 ) 


(2.9) 


Finally, when the liquid fills a slab lying perpendicularly to the longer box edge (L^), the 
free-energy excess becomes: 


APslab = iVi(/i - A) + 2jLl = Py{Vy - v)N + 4^r;2/3/v2/3 ^ 


( 2 . 10 ) 


5 
















whence a pressure of 




( 2 . 11 ) 


In the present analysis, a crossing of free energies as a function of p entails a change in the 
relative stability of two drop shapes. Far from being a hrst-order transition, which can only 
occur in the thermodynamic limit and is not accompanied by any pressure jump, this shape 
transition is an equilibrium hnite-size effect promoted by periodic boundary conditions. In 
fact, all shape transitions become rounded crossovers when thermal fluctuations are taken 
into account (see Refs. 13 E])- Upon equating the free energies two at a time, and providing 
that the sequence of “phases” is the same as found in Ref. [6], we obtain the following 
formulae for the “transition” densities: 



Pcyl—slab Pv U (Pl Pv) • 


( 2 . 12 ) 


In particular, we note that the density range beyond pv where the homogeneous vapor is 


thoroughly stable vanishes in the thermodynamic limit as N This p interval should not 


be confused with the vapor metastability region, which instead is a kinetic concept more 
appropriate to bulk systems. The metastability region of vapor, which extends past the 
T-P coexistence locus inside the liquid region, comprises all nominally-liquid states where 
a quenched vapor system can be maintained gaseous for a long time (that is, much longer 


than the typical observation times) before nucleation of the liquid occurs. 


With a further little effort, one may also derive the size of the liquid droplet in its various 
conformations, obtaining: 



(2.13) 


where the latter quantity represents the thickness of the liquid slab. 

As a hnal comment, we underline that the above theory would also apply with no modifi¬ 
cations to the solid-liquid transition na. However, a testing of the theory against simulation 
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data would actually be impossible in this case, because of the tendency of both phases to 
go deeply metastable (see, e.g., Ref. i), a fact that generally prevents one from observing 
any shape transition. A notable exception is the numerical evidence reported in Fig. 2 of 
Ref. [18], where, thanks to a smart simulation method, the low system dimensionality, and 
very long runs, it was possible to construct a pressure equation of state showing a plateau 
in the density range corresponding to the slab formation. 

III. ASSESSMENT OF THE THEORY: LENNARD-JONES MODEL 

As originally shown in Ref. [6], the sequence of shape transitions in a large periodic cubic 
box is expected to be horn —)■ sph —)■ cyl —)■ slab, then followed by the reversed transition 
sequence with the role of vapor and liquid interchanged. This hnding is conhrmed by the 
results of our simulations for the cut-and-shifted LJ model and gratifyingly reproduced by 
our theory for a = 1 (see Sect.III-A below). The possibility of other stable shapes, which 
as far as we know was never examined before, is discussed in the next Sect. III-B. 

A. Results from sequential simulations 

In a recent paper |1|, we performed extensive Metropolis MC simulations of the cut- 
and-shifted LJ model in the NVT ensemble. This model is characterized by the following 
interaction potential: 



4e [(cr/r)^^ — (n/r)®] — c, for r < rcut 

0 , for r > rcut 


(3.1) 


with Tcut = 2.5(7, where the constant c is chosen in such a way as to make u{r) everywhere 
continuous (from here onward, all quantities will be expressed in the units set by /cB,e, and 
(7, where k-Q is Boltzmann’s constant). The critical temperature of this system is slightly 
less than 1.10 [IH]. In particular, in Ref. |1| we considered the behavior of the system along 
the isotherm T = 0.90, to see whether in the liquid-vapor region equilibrium can be reached 
by traditional simulation methods with local moves. To this end we performed simulation 
runs in a sequence, starting at each density from the last system conhguration produced 
in the previous run at a slightly smaller density. We thus showed that, anywhere within 
the coexistence region, heterogeneous equilibrium can be established in an affordable time. 
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This is proved by the fact that we were able to obtain the known liqnid-vapor coexistence 
densities by integrating the pressnre eqnation of state across the two-phase region. At least 
this was the case for a system of = 1372 particles or smaller, while a larger system of 
4000 particles was fonnd to be plagued by metastability (he., the data collected along the 
forward and backward trajectories were not the same, see Fig. 1 below). 

Having accurate simulation results available in the liquid-vapor region gives the opportu¬ 
nity to make a critical appraisal of the theory presented in Sect. II. For T = 0.90, the liquid- 
vapor coexistence pressure and densities are Py = 0.03146, pv = 0.0451, and p\ = 0.6649, 
respectively |1]. In order to compute we carried out a long NPT simulation of the 
vapor at P = Py, eventually finding an isothermal compressibility of 45.04 (we computed 
Py from the formula = /^((H^) — (H)^)/(H), where /? = and (■■■) is the 

isothermal-isobaric average). The only theory parameter left to set is 7 , and we decided to 
choose it in such a way that Phom-sph roughly coincides with the location of the hrst pressure 
drop in the MC data for N = 4000 [1] (we thus obtained /Sy = 0.19). We report in Fig. 1 
the results from theory, and compare them against MC data. In the top panel, we show 


AF/N — Pv{vv — v) for the homogeneous vapor (cf. Eq. (2.2)) and the various heterogeneous 


phases (cf. Eqs. (2.4), (2.8), and (2.10)). The sequence of stable conformations as a func¬ 
tion of density is as expected, that is horn —sph —)■ cyl —)■ slab. In the middle panel of 
Fig. 1 we compare the theoretical pressure (red line) with MC data from both forward- and 
backward-travelled paths | 1 ]: we see an overall good agreement, especially in the slope of 
the pressure plateaus, but also the location of shape transitions would be well reproduced by 
the theory considering that the data from the backward trajectory are probably the closest 
to equilibrium (indeed, it is much easier for the system to preserve its structure during over¬ 
compression rather than when reducing the density below the transition threshold). Finally, 
the bottom panel in Fig. 1 reports the characteristic size of the liquid drop in each “phase”. 

For high enough densities, other conformations of drop become available to the system, 
as seen in the complete MC pressure equation of state (Fig. 2) where each inflection point 
marks a different shape transition. The theory for these further transitions can be formulated 
following the same steps as before, by everywhere exchanging vapor with liquid, and the 
results are shown in Fig. 3 for the same parameters used to draw Fig. 1. From a glance at 
Fig. 3 we see that the accuracy of the theory is now poorer, though it would be safe to say 
that it still qualitatively reproduces the MC data. A worse agreement with Monte Carlo 






should anyway be expected just for the reason that equilibration is more difficult at high 
density and the spatial dehnition of the interface between liquid and vapor is also poorer. 
In summary, we reached a good agreement between simulation and theory by adjusting one 
single parameter. 

In Ref. |1] we also reported simulation data obtained from sequential runs of = 1500 
Lennard-Jones particles in a periodic cuboidal box with edges in the ratio of 1:1:3, for 
T = 0.90. These results revealed that the spherical “phase” is apparently never stable 
whereas the cylindrical “phase” is strongly reduced in extent. We tested this conclusion by 
the theory of Sect. II and we indeed found that, under compression, the homogeneous vapor 
hrst gives way to a cylindrical drop of liquid, which then changes to a slab upon increasing 
the density further (see Fig. 4). Again, we see a clear correlation between the location of the 
inflection points in the pressure data and the transition thresholds derived from theory. As a 
last note we observe that, according to the same theory, the stability of the cylindrical drop 
progressively reduces upon increasing a, until the cylindrical phase disappears altogether 
(for a ~ 10) and the homogeneous vapor is thereafter directly transformed into the slab 
phase. The use of an elongated box is instrumental to obtaining a wider density interval 
for the study of the planar liquid-vapor interface, which could be useful for example in the 
determination of the surface tension as a function of temperature. 


B. Other drop shapes 

Our next point concerns the evidence, originally found in LJ-model simulations for T = 
0.75 [1], of an additional pressure “plateau” (in fact, a slightly inclined flat region) lying 
between the “cylinder” and the “slab” plateaus. The trick to obtain this result was to start 
the simulation from scratch at every density. We commented in Ref. [1] that, in the interval 
of densities corresponding to the extra plateau, the liquid drop has the shape of a slab with 
a hole. 

We carried out NVT MC simulations of the cut-and-shifted LJ fluid for T = 0.72 in 
a periodic cubic box, for a number of densities close to Pcyi-siab ~ 0.25 (with the values 
of pv and p\ taken from Ref. [20]), always assuming a face-centered-cubic structure for the 
initial conhguration of the run. Clearly, this may not be the best choice to study hetero¬ 
geneous equilibria by simulation, since relaxation times would be much longer than those 
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encountered when performing simulation runs sequentially. However, the point is that by 
this method we may hnd out long-lived metastable states that would not be reached by 
sequential simulations. For each density we hrst equilibrated the system for 5 million cycles, 
then gathering equilibrium statistics for just as many cycles. The pressure data obtained 
for N = 4000 particles were plotted in the middle panel of Fig. 5. In the same hgure we 
also reported the outcome of the theory. We see that the data points are rather clearly 
located over four distinct pressure plateaus, each being representative of a different drop 
conformation of long life. While three of these plateaus have already been identihed as 
representative of spherical, cylindrical, and slab-like drops, in the fourth plateau centered 
at p = 0.25 a visual inspection of the system conhguration reveals that the shape is novel; 
the liquid drop resembles either a double cylinder (DC, see Fig. 6a) or a punched slab (PS, 
see Fig. 6b). Owing to the periodicity of the simulation box, the DC and PS shapes would 
actually be similar: it suffices to move the DC center to a box vertex (and possibly rotate 
the structure just to improve visibility) to realize that a DC is in fact not dissimilar from 
a PS (Fig. 7a). Similarly, by moving the hole center to a box vertex the PS of Fig. 6b ends 
up looking like a DC (Fig. 7b). However, while the boundary of the hole associated with 
a perfect DC is a square, the hole of the PS seen in Fig. 6b is roughly circular; hence, it 
would be wrong to conclude that the two geometries are exactly the same up to a folding 
operation. 

Let us hrst attempt to model the liquid drop as a perfect slab with a cylindrical hole (we 
call this geometric shape “type-1 PS”). Let d and r < be the slab thickness and hole 
radius, respectively. The interface area equals 

A = 2(L^ — Trr^) + 27rrd (3.2) 


with 

V = {Ll- 7rr^)d = Nivi = 17 . 


Upon eliminating r in favor of d, we obtain; 

A{d) = ^ + W^Lld^ - Md). 


The previous equation is correct only provided 0 < r < La;/2, or 


u 

LI 


< d < 


4 U 
4 —71 LI 


(3.3) 


(3.4) 


(3.5) 
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Upon taking x = L^d/Vi — 1 = nr'^d/Vi and y = A/{2L‘l) — 1, the problem is reduced to 
finding the absolute minimum of 

y = -7^ + ^V 7 r(a: + a:2) (3.6) 

in the interval 0 < x < 7r/(4 — tt). The hole only forms if, besides the local minimum 
at X = 0, a deeper (negative) y minimum also occurs (at a certain Xmin corresponding to 
d = dmin oc This requires p to be less than a threshold density (~ 0.175, see Fig. 8 

below; beyond this density, the free energy becomes identical to AFsiab)- With the A so 
determined, the free energy of the PS is given by 


AF = Pv(xv - v)N + 'yA . 


(3.7) 


However, representing the slab hole as a perfect cylinder is too rough an approximation, 
since the physical drop would certainly manage to avoid any sharp edges. A more realistic 
modelization (say, “type-2 PS”) will entail a smooth and curved hole boundary, and the 
most natural solution would be the surface of the innermost half of a torus (IHT, namely 
the part of the torus which lies inside a cylinder having the same symmetry axis as the torus 
and radius equal to the torus major radius - i.e., to the distance r from the center of the 
hole to the center of the tube). While the minor radius of the torus (that is, the radius 
of the tube section) has to be half of the slab thickness (d/2), the major radius must obey 
d/2 < r < Lx/2 (when r = d/2 the hole closes and the torus becomes a horn torus). 

We used Pappus’ centroid theorem to compute the area and volume of the IHT of radii 
P< = d/2 and i?> = r. From the general formulas, 

4 

Aiht = 27r^P<P> — 47rP| and VjuT = , (3.8) 

o 


it follows that the A and V in Eqs. (3.2) and (3.3) should be replaced with 

1 


A = 2(L^ — nr^) + n'^dr — vrd^ and V = {Li — 7rr^)d -|- ^Tr^d^r — ^d^ . (3.9) 

4 6 


Using the equation V = V\ to simplify the A expression, we eventually obtain; 


A 







(3.10) 


where r is the solution to 

2 1 Vl 1^2 n 

F - -7idr H- - + -d^ -^ = 0 . 

4 vrd 6 TT 


(3.11) 
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Should the solutions of Eq. (3.11) be both valid, the one must be chosen that provides the 


minimum A value. Then, the free energy follows from the general Eq. (3.7). 


Finally, we considered a drop having the shape of two equal cylinders crossing each other 
at right angles. We call this shape a DC (also the cylinder axes must intersect with one 
another). While the length of the axes is the cylinder radius r should be consistent 

with the V = Vi condition. In order to compute the area and volume of a DC, it comes 
useful to consider the solid body which the two cylinders have in common, which goes under 
the name of hicylinder (or Steinmetz solid). Its area and volume can easily be obtained by 
multiple integration, and the results are 16r^ and (16/3)r^. Clearly, the area and volume of 
a DC then read: 


16 

A = 4:7irLx — 16r^ and V = 27ir‘^Lx — —r^ 


(3.12) 


From V = Vi a. third-order algebraic equation is obtained for r, which turned out to have 
only one solution satisfying 0 < r < Upon plugging this r in A, we again derive AF 


from Eq. (3.7). 


For T = 0.72, the free-energy curves for the three just considered shapes were plotted as 
a function of the density in the top panel of Fig. 8. In the same picture, the lower envelope 
of the free energies for the common shapes was reported for comparison (red line). We see 
that none of the PS nor the DC ever provide the shape of the stable drop. However, a type-2 
PS with r = Lxl2 is nearly stable close to Pcyi-siab — 0.25, which may explain why this drop 
shape was observed in our simulations from scratch (as seen in Fig. 8, suitable type-1 PS and 
DC also exist that offer not too bad solutions near Pcyi-siab)- Beyond p = 0.24, the optimal 
type-2 PS has r = d/2, meaning that the hole eventually closed and the torus became a 
horn torus (however, this degenerate type-2 PS is obviously no longer reliable as a drop 
shape). Finally, notice that changing T = 0.72 to T = 0.90 does not modify the theoretical 
picture in any respect. The conclusion is that, according to the present theory, the three 
novel shapes considered in this section are metastable. Hence, the sequence of stable drops 
for T = 0.72 remains the same identihed in Sect. HI-A for T = 0.90 and a = 1. 

IV. CONCLUSIONS 


Below the critical temperature, the structure of a heterogeneous fluid simulated under pe¬ 
riodic boundary conditions is sensitive to the imposed density. Each change of conformation 
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observed within the binodal line goes along with a characteristic fall in the pressure. In this 
paper, the by-now classical evidence of spherical, cylindrical, and slab-like liquid drops was 
re-examined from the viewpoint of a theory having its roots in a proof-of-concept calculation 
found in Ref. [5|. In spite of its simplicity, this theory well reproduced the behavior of a 
Lennard-Jones fluid along the liquid-vapor region, for both a cubic and an elongated box. 

Further drop shapes emerged near the crossover region from cylindrical to slab-like when 
simulations, rather than being performed in a sequence as is common practice, were started 
for each density from a face-centered-cubic conhguration. In this case, the liquid drop 
occasionally exhibited the shape of a punched slab with a roughly circular hole. Employing 
our theory, we found that this peculiar shape of drop is actually only metastable, he., it is 
a long-lived conformation which however is doomed to decay. 
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FIG. 1. (Color online) Comparison between theory and MC data [1| from sequential simulations 
of the cut-and-shifted LJ model in a periodic cubic box (A^ = 4000 and T = 0.90). In order to 
equilibrate the system 10® MC cycles were produced at each state point, followed by other 4 x 10® 
cycles over which equilibrium averages were computed. Top: AF/N — Pv{vv — v) for the various 
competing “phases” (green, “horn”; cyan, “sph”; blue, “cyl”; black, “slab”). The lower envelope of 
all the free-energy curves (thick red line) marks the equilibrium curve. The vertical lines mark the 


location of the shape transitions (see Eqs. (2.12)). Middle: system pressure relative to P^- The dots 
are the MC data from Ref. [3] (open dots, forward trajectory; full dots, backward trajectory) while 
the thick red line is the equilibrium pressure according to our theory. Bottom: the equilibrium 


drop size (thick red line, see Eqs. (2.13)) vs. half box size (black line). 
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FIG. 2. MC pressure data for the cut-and-shifted LJ model in a periodic cubic box with N = 4000 
and T = 0.90 (only the forward trajectory is shown). The horizontal line lies at the level of Pv 
The vertical lines mark the liquid-vapor coexistence densities, and p\. 
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FIG. 3. (Color online) Comparison between theory and MC data from sequential simulations of 
the cut-and-shifted LJ model in a periodic cubic box [N = 4000 and T = 0.90). Same as Fig. 1 
but the high-density region is shown here. 
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FIG. 4. (Color online) Comparison between theory for a = 3 and MC data from sequential 
simulations of the cut-and-shifted LJ model in a periodic cuboidal box (N = 1500, T = 0.90, and 
= 0.175). In strike contrast with the a = 1 case, the spherical drop would never be stable for 
a = 3. The cylindrical drop is only observed in a narrow range of densities close to 0.1 whereas 
the slab-like drop is stable over a much wider density interval. 
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FIG. 5. (Color online) MC pressure data for the cut-and-shifted LJ model in a periodic cubic box 
[N = 4000 and T = 0.72) vs. theory (for this case, = 0.0062, = 0.0086, p\ = 0.7722, j3^ = 

0.78 |20j . whereas the value of 173.62, resulted from a long NPT simulation of the liquid at 
Pv)- The data points (full dots in the middle panel) were obtained independently from one another, 
by starting the simulation from scratch at each density. We see that each point lies on a theoretical 
line, except for the points around p = 0.25 which appear to lie on a different plateau (dashed line). 
In the corresponding range of densities the liquid drop resembles either a DC or a PS. 
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FIG. 6. (Color online) Cut-and-shifted LJ model for N = 4000 and T = 0.72; results from MC 
simulations started from scratch, (a) A system snapshot taken at the density p = 0.24; (b) another 
snapshot taken at p = 0.26. In these pictures, each particle was given an effective diameter a. The 
symbol colors were chosen according to the number of nearest neighbors each particle has, in turn 
defined as the number of particles within a distance of 1.45 ct from the given particle (1.45(T being 
the position of the first non-zero minimum of the radial distribution function of the liquid at pi). 
While the liquid drop in (a) resembles a DC, the drop in (b) looks like a PS. 
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FIG. 7. (Color online) Same as in Fig. 6, but after the particle coordinates have been transformed 
as explained in the main text. The outcome was that a DC became a PS, and vice versa. 
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FIG. 8. (Color online) Shape transitions for the cut-and-shifted LJ model in a periodic cubic box 
{N = 4000 and T = 0.72): theoretical results (see parameters in Fig. 5 caption). Top: AF/N — 
Pv{vy — v) for various competing “phases” (blue full line, “cyl”; black full line, “slab”; green long- 
dashed line, type-1 PS; cyan long-dashed line, DC; magenta long-dashed line, type-2 PS - notice 
that the free-energy curves for “horn” and “sph” are not shown). The thick red line marks the 
equilibrium curve (that is, the lower envelope of all the free-energy curves, including “horn” and 
“sph”). The vertical lines mark the location of the shape transitions. Bottom: the equilibrium 
drop size (thick red line) vs. half box size (black full line). The dotted and dashed lines respectively 
refer to the optimal d and r values for the two types of PS and the DC. 
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